%% produces simulations of firms over lifcycle for Table 6
m_minus = zeros(11,11);
m_end   = zeros(11,11);
b_end   = zeros(11,11);
y_end   = zeros(11,11);
n_end   = zeros(11,11);

for year=1:11

    T = (11-year+1)/dt;%number of periods to simulate the firm
    year_age0 = (year-1)/dt;
    n_sim = ceil(n_sim_base*entry_scaled(year));

    b_hist =zeros(n_sim,n_years/dt);
    va_hist=zeros(n_sim,n_years/dt);
    n_hist =zeros(n_sim,n_years/dt);
    i_alive=zeros(n_sim,n_years/dt);
    i_exit =zeros(n_sim,n_years/dt);

    for ii=1:n_sim
        bd=rand;
        z0=zlow*(bd<=q_hig)+zmid*(bd>q_hig)*(bd<=q_hig+q_mid)+zhig*(bd>q_hig+q_mid);
        inc_draw=rand;
        %draw leverage of firms at birth before and after the subsidy
        if year<8
            b0=b_up_ss*(bd<=q_hig)+b_md_ss*(bd>q_hig)*(bd<=q_hig+q_mid)+b_dn_ss*(bd>q_hig+q_mid);
        elseif year ==8
            if inc_draw>incid/2
                b0=b_up_unfund*(bd<=q_hig)+b_md_unfund*(bd>q_hig)*(bd<=q_hig+q_mid)+b_dn_unfund*(bd>q_hig+q_mid);
            else
                b0=b_up_fund*(bd<=q_hig)+b_md_fund*(bd>q_hig)*(bd<=q_hig+q_mid)+b_dn_fund*(bd>q_hig+q_mid);
            end
        else
            if inc_draw>incid
                b0=b_up_unfund*(bd<=q_hig)+b_md_unfund*(bd>q_hig)*(bd<=q_hig+q_mid)+b_dn_unfund*(bd>q_hig+q_mid);
            else
                b0=b_up_fund*(bd<=q_hig)+b_md_fund*(bd>q_hig)*(bd<=q_hig+q_mid)+b_dn_fund*(bd>q_hig+q_mid);
            end
        end

        t=year-1;inde_alive=1;t_sim=0;age=0;
        while t_sim<T && inde_alive==1
            t=t+dt;t_sim=t_sim+1; age=age+dt;

            if t<8.5
                b_bar=b_bar_ss;va=A_ss;
            else
                b_bar=b_bar_iss;va=A_iss;
            end
            if year<8 && abs(t-8.5)<0.01,b0=b0*A_ss/A_iss;end

            if b0<b_bar
                ell_t       =  fnval(cub_spl_l,[max(0,b0);max(va,Amin)]);
                drift_I_t   =  fnval(cub_spl_d,[max(0,b0);max(va,Amin)]);
                drift_dlb   = dt*(ell_t/b0+sigma^2/2-rho-drift_I_t);
                drift_dlz   = dt*(drift_I_t-sigma^2/2);
                dw  = randn;
                l_b = log(b0) + dt*drift_dlb-sqrt(dt)*sigma*dw;
                b   = exp(l_b);
                l_z = log(z0) + dt*drift_dlz+sqrt(dt)*sigma*dw;
                z   = exp(l_z);
                b0=b;z0=z;
                b_hist(ii,year_age0+t_sim) =b0*z0*va;
                va_hist(ii,year_age0+t_sim)=va*z0;
                n_hist(ii,year_age0+t_sim) =va*z0*(eps/(eps-1)*w_t(year_age0+t_sim))^(-eps);
                i_alive(ii,year_age0+t_sim)=1;
            else
                draw_neg=rand;
                if  draw_neg>phi
                    inde_alive= 0;
                    i_exit(ii,year_age0+t_sim:end) =1;
                else
                    b0=b_bar*a;
                    b_hist(ii,year_age0+t_sim) =b0*z0*va;
                    va_hist(ii,year_age0+t_sim)=va*z0;
                    n_hist(ii,year_age0+t_sim) =va*z0*(eps/(eps-1)*w_t(year_age0+t_sim))^(-eps);
                    i_alive(ii,year_age0+t_sim)= 1;
                end
            end
        end
    end

    %take sum year
    for age=1:11-(year-1)%age
        t=year-1;%initial year
        m_minus(age,t+age) = m_minus(age,t+age)+sum(i_alive(:,year_age0+(age-1)/dt+1));
        m_end(age,t+age)   = m_end(age,t+age) +sum(i_alive(:,year_age0+age/dt));
        b_end(age,t+age)   = b_end(age,t+age)+sum( b_hist(:, year_age0+age/dt));
        y_end(age,t+age)   = y_end(age,t+age)+sum(va_hist(:, year_age0+age/dt));
        n_end(age,t+age)   = n_end(age,t+age)+sum( n_hist(:, year_age0+age/dt));
    end
end

%% stack
age =0:10;
year=2010:2020;

age_mat= age'*ones(1,11);
year_mat=ones(11,1)*year;

m_minus_vec = reshape(m_minus,11*11,1);
m_end_vec   = reshape(m_end,11*11,1);
b_end_vec   = reshape(b_end,11*11,1);
y_end_vec   = reshape(y_end,11*11,1);
n_end_vec   = reshape(n_end,11*11,1);
age_end_vec = reshape(age_mat,11*11,1);
year_end_vec= reshape(year_mat,11*11,1);

m_minus_final=m_minus_vec(m_minus_vec>0);
m_end_final  =m_end_vec(m_end_vec>0);
b_end_final  =b_end_vec(b_end_vec>0);
y_end_final  =y_end_vec(y_end_vec>0);
n_end_final  =n_end_vec(m_end_vec>0);
age_end_final=age_end_vec(m_end_vec>0);
year_end_final=year_end_vec(m_end_vec>0);

data_mat = [m_minus_final,m_end_final,b_end_final,y_end_final,n_end_final,age_end_final,year_end_final];